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We adapt a previous model and analysis method (the master stability function), extensively used 
for studying the stability of the synchronous state of networks of identical chaotic oscillators, to 
the case of oscillators that are similar but not exactly identical. We find that bubbling induced 
desynchronization bursts occur for some parameter values. These bursts have spatial patterns, which 
can be predicted from the network connectivity matrix and the unstable periodic orbits embedded 
in the attractor. We test the analysis of bursts by comparison with numerical experiments. In the 
case that no bursting occurs, we discuss the deviations from the exactly synchronous state caused 
by the mismatch between oscillators. 
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I. INTRODUCTION 
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In this paper we study the synchronization of networks 
of coupled chaotic units that are nearly, but not exactly, 
identical. In particular, we will be concerned with the 
spatial patterns of desynchronization bursts that appear 
when this synchronization is present but intermittent. 

When two or more identical dynamical systems are 
coupled, they can synchronize under appropriate circum- 
stances. The synchronization of chaotic units has been 
studied extensively Q, 01_and is of significance in biol- 
ogy laser physics 0-0, an d other areas [IollTT|. 
At the same time, the importance of complex networks 
has been recently appreciated, and progress has been 
made towards their understanding, including character- 
istics that mi ght help distinguish qualitatively different 
networks [T^-|l4). The dynamics of a network of cou- 
pled oscillators, and, in particular, its synchronization, 
has therefore emerged as a subject of great interest. 

Pecora and Carroll have proposed a model and 
analysis method (the master stability function) for the 
study of the stability of the synchronous state of networks 
of identical coupled chaotic units, and this technique has 
recently been extensively applied 0, to study the 
synchronization properties of different kinds of networks 
of identical noiseless chaotic units. These networks in- 
clude small world 0] and scale- free networks . 

The analysis of network synchronization by use of the 
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master stability function technique has so far assumed all 
the units to be identical and noise-free, so that an exact 
synchronized state is possible. In practice, however, even 
if one strives to make the oscillators the same, they are 
still expected to have a small amount of parameter mis- 
match, and a small amount of noise is also expected to be 
present. Under such circumstances, it is known that the 
synchronization can be interrupted by sporadic periods 
of desynchronization (bursts). The bursts are typically 
caused by a periodic orbit that is embedded in the syn- 
chronized chaotic attractor and is unstable in a direction 
transverse to the synchronization manifold. This phe- 
nomenon is commonly referred to as bubbling [l9 | -|2l | . 
and has been studied extensively for two coupled oscilla- 
tors pll^. 

Our purpose in this paper is to study desynchroniza- 
tion bursts in networks of coupled chaotic nonidentical 
units. (Noise has a similar effect but will not be treated 
in this paper.) We will use the master stability function 
approach and, in order to account for the possibility of 
bubbling, we will also extend this approach to include 
the stability of embedded periodic orbits. In this case, 
the bursts have the added feature of having spatial pat- 
terns on the network, and we find that these patterns 
can be predicted from the network connectivity matrix. 
We will show how these bursts affect different parts of 
the network in different ways. In particular, we will see 
how adding connections in a ring can destabilize precisely 
those nodes that are the most connected, leaving other 
parts of the network substantially synchronized. (This a 
somewhat counterintuitive effect related to the fact that, 
in some cases, increasing the coupling strength destabi- 
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lizes the synchronous state |l5|.[2*4|.) 

Arbitrarily small amounts of mismatch will eventually, 
through the bubbling mechanism, induce desynchroniza- 
tion bursts. We will show that some of the spatial pat- 
terns of this possibly microscopic mismatch might get 
amplified to a macroscopic size in the bursts. We will 
discuss how one can use knowledge of the parameter mis- 
match of the dynamical units in the network to decrease 
the effective size of the mismatch driving the bursts, 
thereby improving the robustness of the synchronization. 

If synchronization is desired, the network and the pa- 
rameters should be constructed so that the synchronous 
state for the identical oscillator system is robustly stable 
(this implies the absence of noise or mismatch induced 
desynchronization bursts). Even then, the synchroniza- 
tion will not be perfect if the oscillators have parameter 
mismatch. We will describe the characteristics of the de- 
viations from exact synchronization in terms of the mis- 
match and the master stability function. 

This paper is organized as follows. In Sec. II we re- 
view the master stability function approach and apply 
it to the case of coupled Rossler units. We also discuss 
the bubbling mechanism by including the embedded pe- 
riodic orbits in the master stability function analysis. In 
Sec. Ill we numerically consider particular networks as 
examples and show the resulting bursts and their spa- 
tial patterns. The patterns we obtain are long and short 
wavelength modes in a ring and localized bursts produced 
by strengthening of a single connection in a ring. In 
Sec. IV we study the effects of the spatial patterns of the 
mismatch in the development of the bursts. In Sec. V we 
study the deviations from the synchronous state caused 
by the mismatch when the synchronous state of the iden- 
tical oscillator system is stable. In Sec. VI we summarize 
our conclusions. 



II. MASTER STABILITY FUNCTION AND 
BUBBLING 

We now briefly review the master stability function 
approach introduced in . Consider a system of N dy- 
namical units, each one of which, when isolated, satisfies 
Xi = F(Xi,fii), where i = 1,2,.. . N, and Xi is the d- 
dimensional state vector for unit i. In [la ] the parameter 
vectors fa are taken to be the same, fa = ji. Here, how- 
ever, the parameter vectors fa are in general different for 
each unit, but we assume the difference, or mismatch, 
between them to be small. Generalizing the situation 
treated in Ref. 15] to the case where the individual units 
are not identical (i.e., the fa are not all equal), the system 
of coupled dynamical units is taken to be of the form 

N 

X t = F(X U fa) -gJ2 GijH(Xj), (1) 
i=i 

where the coupling function H is independent of i and j, 
and the matrix G is a Laplacian matrix . Gij = 0) de- 



scribing the topology of network connections. For i =/= j, 
the entry Gy is zero if oscillator i is not connected to os- 
cillator j and nonzero otherwise. The nondiagonal entries 
of G are determined by the connections, and the diagonal 
elements are the negative of the sum of the nondiagonal 
matrix elements in their row. The coupling constant g 
determines the global strength of the coupling. 

Assume first that all the dynamical units are identical, 
that is, fa = fi. We will refer to this situation as the ide- 
alized case. In this case there is an exactly synchronized 
solution Xi = X 2 = ■ ■ ■ = Xn = s(t) whose time evo- 
lution is the same as the uncoupled dynamics of a single 
unit, s — F(s,/i). This convenient result arises because 
the Pecora-Carroll model uses the particular choice of 
coupling in JQ) that ensures that the summation is iden- 
tically zero when all of the Xj are equal. We will denote 
this synchronization manifold, X\ = X 2 = • ■ ■ = Xn, by 
M. This manifold is a d - dimensional surface within the 
Nd - dimensional phase space of Eq. • 

The stability of the synchronized state can be deter- 
mined from the variational equations obtained by con- 
sidering an infinitesimal perturbation q from the syn- 
chronous state, Xi(t) — s(t) + €i(t), 

N 

ei = DF(s)e t - g^G ij DH{s)e j . (2) 
i=i 

Let e = [ei, £2, . . . , e_/v] be the d x N matrix representing 
the deviation of the entire network from the synchronized 
state. In matrix notation, Eq. (|2J) becomes 

e = DF{s)e- gDH{s)eG T . (3) 

While © allows for nonsymmetric coupling, we hence- 
forth assume the coupling matrix G to be symmetric, 
G = G T . We write the symmetric matrix G as G = 
LAL T , where A is the diagonal matrix of real eigenval- 
ues Ai, A2, . . . , Ajy of G and L is the orthogonal matrix 
whose columns are the corresponding real orthonormal 
eigenvectors of G (L T L — I). Define the d x N matrix 
V = fail T2j • • • j Vn] by e = r/L T . Then Eq. © is equiva- 
lent to 

7) = DF(s)rj - gDH(s)T]A. (4) 
Componentwise, 

r )k = {DF{s)-g\ k DH(s))r ]k . (5) 

The quantity rjk is the weight of the k th eigenvector of G 
in the perturbation e. The linear stability of each 'spatial' 
mode k is determined by the stability of the zero solution 
of JHJ). As a consequence of the condition Y2j Gij = 0, 
there is a special eigenvalue, A = 0, whose eigenvector 
is ejsi — [1,1,1,...,!], corresponding to perturbations 
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in the synchronization manifold M. Since these are not 
perturbations from the synchronous state, the analysis is 
focused on the perturbations corresponding to nonzero 
eigenvalues. 

By introducing a scalar variable a = gXk, the set of 
equations given by JSJ) can be encapsulated in the single 
equation, 

f) = (DF(s)-aDH(s))r). (6) 

The master stability function *f>(a) [l5| is the largest Lya- 
punov exponent for this equation for a typical trajectory 
in the attractor. This function depends only on the cou- 
pling function H and the chaotic dynamics of an individ- 
ual uncoupled element, but not on the network connectiv- 
ity. The network connectivity determines the eigenvalues 
Afe (independent of details of the dynamics of the chaotic 
units). In the sense of typical Lyapunov exponents, the 
stability of the synchronized state of the network is de- 
termined by = sup fc ty(g\k), where ^F* > indicates 
instability. Thus the Pecora-Carroll model cleanly breaks 
the stability problem into two components, one from the 
dynamics (obtaining 'J (a)) and one from the network 
(determining the eigenvalues A&). 

In contrast to previous work using the master stabil- 
ity function technique, in this paper we are interested 
in the dynamics of systems in which a small parameter 
mismatch is present. (Even though in this paper our 
examples are restricted to the case of mismatch, we em- 
phasize that the same type of bursting phenomenon is 
expected for identical oscillators if noise is present [T^| - 
|22|.') Although the synchronization manifold M present 
in the dynamics of the idealized system is, in general, not 
invariant for the system with mismatch, it still may pro- 
vide a useful approximation to the dynamics in systems 
with small mismatch. If M is stable for the idealized 
system, and the mismatch is small enough, then trajec- 
tories near M will tend to stay near M , and we regard the 
vicinity of M to be the "synchronized" state. However, 
stability of M in the idealized case of identical oscillators 
is not sufficient to guarantee robust synchronization in 
a real system where the oscillators are not identical [l9j- 
[22j| . While in the vicinity of the synchronization man- 
ifold M, a typical trajectory will eventually follow very 
closely a periodic orbit embedded in the attractor of the 
idealized system. Some of these periodic orbits may be 
unstable in a direction transverse to M. When in the 
vicinity of a transversally unstable periodic orbit, mis- 
match (or noise) will cause the trajectory have a compo- 
nent in the direction transverse to M and hence to leave 
the vicinity of the synchronization manifold M. If there 
are no other attractors, the trajectory will eventually re- 
turn to the vicinity of M, and the process will repeat, 
the result being bursts of desynchronization sporadically 
interrupting long intervals of near synchronization. This 
type of dymanics is called bubbling [lj| . 

Thus, in the presence of mismatch (or noise), to deter- 
mine the robustness of synchronization, it is necessary to 



determine the transverse stability of the embedded pe- 
riodic orbits for the noiseless system of identical oscilla- 
tors. For coupling as in J[J, this analysis is independent 
of the network, and such analyses have been carried out 
before, e.g., for the analysis of two coupled oscillators in 
Ref. [23j ■ Equation © can be used as before to construct 
the master stability function for each periodic orbit, if the 
appropriate periodic trajectories are inserted for s(t) in 
©. 

As an example, in this paper we work with the Rossler 
system |2q : 

x = -(y + z), (7) 
y = x + ay, 
z = b + z{x — c). 

In terms of our previous notation, d = 3, fj, = [a,b,c] T , 
and X = [x,y, z] T . We choose the parameters of the 
idealized system to be a = b = 0.2, c = 7. For these pa- 
rameters, the system has a chaotic attractor (see Fig.^l. 
We found the periodic orbits embedded in this attractor 
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FIG. 1: Rossler attractor (projection onto x — y plane) and 
embedded period 1 orbit, displayed as a thick white curve 
inside the attractor. The parameters are a — b — 0.2, c = 7. 



up to period five, and performed the analysis described 
above on them. We found these orbits by looking at the 
Poincare surface of section {y = 0,ir<0}. To a good 
approximation, in this surface of section the dynamics is 
well described by a one dimensional map a; n +i = f(x n ), 
which we approximated using a polynomial fit. From 
this approximation to /, we determined periodic orbits 
of period p by using Newton's method to find the roots of 
x = f p (x), where f p denotes the p times composition of 
/. We found one period 1 orbit, one period 2 orbit, two 
period 3 orbits, three period 4 orbits, and four period 5 
orbits. Using coupling through the x coordinate, 

H([x,y lZ ] T ) = M,0] T , (8) 
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we obtained a stability function '5 (a) for each orbit, the 
largest of which will determine if the synchronization is 
robust. Results are shown in Fig.0 For all values of a, 
we found that the master stability function corresponding 
to the period 1 orbit (thick dashed curve) is larger than 
that for a typical chaotic orbit (thick continuous curve), 
as well as those for the other periodic orbits we have 
found (several of which are shown as thin curves). 



III. EXAMPLES 

In this Section we provide examples of spatially pat- 
terned bursting by considering different configurations of 
the chaotic units. We will first work with the units con- 
nected in a ring with each connection of equal strength. 
The Laplacian matrix G for this arrangement is 
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FIG. 2: Master stability function 9(a) for atypical trajectory 
in the attractor (thick continuous curve), for the period 1 orbit 
(thick dashed curve), and for periodic orbits up to period 4 
(thin curves). The curves for the four period 5 orbits are 
similar to the latter and were left out for clarity. 



G = 



/ 2 -1 ■■• -1 
-1 2 -I ■■• 
-I 2 -I ••• 



V-i o 



(9) 
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and its eigenvalues are given by Xk = 4 sin (^) 
A/c = Ajv-is, each eigenvalue has multiplicity two, with 
the exception of A at = 0, and, if N is even, An = 4. 
The matrix G is shift invariant, that is, its entries satisfy, 
modulo N , Gij — Go,i-j- Under these conditions, the di- 
agonalization procedure described above corresponds to a 
discrete Fourier transform j24[. For the eigenvalue we 
choose the eigenvector Wk given by Wk oc [sin(^-)]j^ 1 

for 1 < k < f , and by w k oc [cos(^)]f =1 for f < 
k < N. (Due to the degeneracy of the eigenvalues in this 
case, there is some arbitrariness in choosing the eigen- 
vectors.) Thus, the longest wavelength modes have the 
smallest eigenvalues, and viceversa. 



A. Long wavelength burst 



Based on the discussion above, bubbling induced burst- 
ing should occur whenever the master stability function 
for a typical chaotic orbit in the attractor is negative 
for a — g\k and all k, while the period one orbit has 
positive master stability function for a — gXu for some 
value of k. Denoting the master stability function for a 
typical chaotic orbit by 9o(a) (Thick continuous curve 
in Fig. an d for the period one orbit by 'J'i(a) (Thick 
dashed curve in Fig. , the bubbling region of a corre- 
sponds to Wo (a) < 0, 'i'i(a) > 0. In our example, this 
region corresponds to 0.16 < a < 0.48 or 3.8 < a < 4.5. 
The range 0.48 < a < 3.8 will be referred to as the stable 
region, and the remaining zone will be called the unstable 
region. 

If a network of slightly mismatched chaotic systems 
coupled according to ^s to be robustly synchronizable 
without bursts of desynchronization, gX k must lie in the 
stable region for all k, where Xk is the fcth eigenvalue of 
G. If gXk lies in the stable region for some k and in the 
bubbling region for other k, then bubbling will typically 
occur. 



First we consider a case in which bursting of the longest 
wavelength mode occurs. We consider N = 12 and g = 
0.71. With these values, the longest wavelength mode 
corresponds to a = gX\ ~ 0.19. This value is in the 
bubbling region, and all other modes are in the stable 
region. 

To introduce heterogeneity in the dynamical units, we 
imagine that we have mismatch predominantly in one of 
the parameters, say a. We simulate this mismatch by 
adding random perturbations to the parameter a of each 
oscillator. These perturbations are uniformly distributed 
within a ±0.5% range; i.e., ai is chosen randomly in the 
interval [0.995a, 1.005a], where a is the parameter value 
of the unperturbed system (a = 0.2). The parameters 
b and c were taken to be the same for each oscillator, 
bi = b = 0.2, Ci = c = 7. How a particular choice of the 
mismatch affects the bubbling process will be discussed 
in Section IV. 

We solved the 12 coupled differential equations 
(Eq. (Q) with the initial conditions chosen near the at- 
tractor in the synchronization manifold. In Fig. [3] we 
plot the quantity xi - x 6 for 1000 < t < 1600. Most 
of the time, this variable is close to zero, as expected 
if the oscillators are synchronized. Approximately at the 
time t = 1380, this difference grows, reaching magnitudes 
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FIG. 3: xi - 16 as a function of time for N = 12 Rossler 
systems connected in a ring with g = 0.71. Note the desyn- 
chronization burst which starts at t « 1380. 



close to 3. By time t = 1500, the difference has decreased 
and is again close to zero. 



The desynchronization burst can be observed developing 
mainly at the longest possible wavelength. 
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FIG. 5: Xj — Xj-i versus the node index j for t = 1360 
(open triangles), t = 1385 (open circles), and t = 1410 (open 
squares). Note that the burst is absent first and grows with 
a long wavelength pattern. 




FIG. 4: xi versus yi for 1372 < t < 1392. During this period, 
which corresponds approximately to the starting point of the 
burst in Fig. [3] the trajectory follows closely the transversally 
unstable period 1 orbit embedded in the attractor (See Fig.0. 



To confirm the mediating role of the embedded unsta- 
ble periodic orbits in the development of the desynchro- 
nization burst, we show in Fig. ^a plot of x\ versus y\ 
from t = 1372 to t — 1392, which is near the start of the 
burst. During this time, the trajectory closely follows the 
period 1 orbit, which is the most transversally unstable 
of the periodic orbits. Similar observations have been 
previously reported for two coupled chaotic systems |23j . 

Finally, in Fig. [S] we plot X j " Xj i cLS cL function of 

j, the oscillator index, for t = 1360 (open triangles), 
t = 1385 (open circles), and t — 1410 (open squares). 



When subsequent bursts were studied in the same way, 
it was found that the phase of the long wavelength burst 
assumed only one value. This is due to the fact that the 
mismatch is 'frozen', that is, each oscillator has a given 
set of parameters which differs by a given amount from 
the mean values. This fixed spatial heterogeneity favors 
certain spatial patterns over others. We will discuss this 
in more detail in Section IV. 



B. Short wavelength burst 

Short wavelength bursting can be expected, for exam- 
ple, when N = 8 and g = 1.09. In this case the value of 
Afe corresponding to the shortest wavelength mode yields 
g\k = 4.36, which is in the bubbling region, while all 
the other modes are in the stable region. In this case the 
observation of the bursts is more difficult, as the transver- 
sal instability of the orbits and the transversal stability 
of the attractor are less pronounced (compare ^(4.36) for 
this case vs. ^(0.19) for the previous example in Fig.EJ- 
Accordingly, the perturbations of the parameter a were 
made larger, with perturbations randomly chosen with 
uniform density within a ±6% range of the ideal values 
of the parameter (a = 0.2). In principle this is not nec- 
essary, as a burst will eventually occur after long enough 
time. In practice, however, it is necessary to reduce the 
waiting time to a reasonable value. As before, the cou- 
pled equations were solved with an initial condition on 
the synchronization manifold. In Fig. we show y 2 — y\ 
as a function of time for one choice of initial conditions. 
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FIG. 6: 2/2 — yi as a function of time for 8 Rossler systems in 

a ring. The coupling strength g was 1.09. The desynchroniza- FIG. 7: Vj-yj-i versus the node index j for t = 17970 (open 

tion burst develops at t » 18000, although it is not as sharp triangles), t = 18270 (open circles), and t = 18570 (open 

due to the smaller magnitude of the transversal Lyapunov squares). The desynchronization burst has a short wavelength 

exponents (*(4.36) in Fig.0. spatial pattern. 



The difference y2 — yi is usually positive and of magni- 
tude close to 1. This asymmetry is not a surprise since 
the oscillators are slightly different. For the relatively 
large value of the mismatch used, this is the "synchro- 
nized state". It is seen in Fig^that the difference y^ — yi 
increases rapidly at around t w 18400, and soon reaches 
values higher than 7. It remains large for a longer time 
than in the case of the long wavelength burst (see Fig. EI) 
and decays more slowly as well. This is in qualitative 
agreement with the smaller absolute values of the master 
stability functions for the short wavelength mode, both 
for typical orbits on the attractor and for the periodic 
orbits. 

In Fig. [7] we plot yj — yj-i as a function of j, the os- 
cillator index, for t = 17970, t = 18270 and t = 18570. 
As expected, the burst mainly affects the shortest wave- 
length mode. In order to perform a more quantitative as- 
sessment of these observations, we define = {( 



where 



([m-k]x) 2 }i for 1 < k < f and £f 

[Vk]x is the x component of the three dimensional vector 
r/k- The quantity is proportional to the magnitude of 

the vector \r)\j, sin (—fj^J +i]i,N-k cos ^ ne com- 

ponent in the Fourier decomposition of the perturbation 
e associated with modes k and N — k. Thus, the quan- 
tity £fc represents the weigth of the modes associated to 



the eigenvalue which have, for 1 < k < ^, a wave- 
length of -g'th of a full wavelength. In Fig. we plot as 
a function of time the quantities for k — 1, 2, 3, 4. The 
short wavelength mode (k — 4, upper curve) is dominant 
during the burst. 




18000 18500 19000 19500 

t 



FIG. 8: Cfc as a function of time for k = 1,2,3,4. The 

shortest wavelength component corresponds to k — 4 (top 

curve). The second shortest one is k = 3 (middle curve). The 
curves corresponding to k = 1, 2 are close to zero. 



C. Localized burst 

In the above examples all links had equal weights. As 
an example of a case with unequal link weights we con- 
sider the case where the previous network is modified 
by doubling the strength of one of the links. Let the 
link whose strength is doubled be the link that connects 
nodes p and p + 1. For example, for p = 4, N — 8, this 
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yields the Laplacian matrix 
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(10) 



Adopting the analysis technique of Ref. |2g, we can 
show that such an enhanced connection has the conse- 
quence that the largest eigenvalue of G corresponds to 
an eigenfunction that is exponentially localized to the 
region near the strong connection. That is, for large N, 
the components of this eigenfunction decay exponentially 
as the distance between the localized region and the node 
corresponding to a component increases. Using the ideas 
of Ref. 26], we now provide this analysis. The equations 
for the eigenvector w and eigenvalue A are 



2w r , 



P-1 

f 3w v 



+ 3w p = \w. 



u p+ 2 - 2w p + 3w p+ i = Xwp+i, 
— Wj-i — Wj+x + 2u>j = XtVj, 



(11) 



for, respectively, nodes p, p + 1 and j different from p 
or p + 1. Since nodes p and p + 1 are identical, we can 
assume, for k > 0, w p+ i + k = ±Wp-k- Furthermore, we 
propose the ansatz w p+ i + k oc t , for k > 0. This will 
be a good approximation if the mode is localized (i.e., 
\t\ < 1), and the network is big enough that <C 1. 
In the antisymmetric case, u> p +i+fc = — u> p _fc, Eqs. 
yield. 



t = A, 
2 = A 



(12) 



which gives 



-t - tr L + 



t- 1 A - 16 
~3 ' ~ T' 



(13) 



Compare this eigenvalue with the largest eigenvalue for 
the network in which all links have equal strength, which 
has a value of 4. The symmetric solution, w p+ i + k = 
Wp-k, yields t — 1 and A = 0, corresponding to the eigen- 
vector [1, 1, ... 1] of perturbations in the synchronization 
manifold. The smallest nonzero eigenvalue remains un- 
changed. 

As an example, we show the localized desynchroniza- 
tion bursts produced by one of these strengthened con- 
nections for the case N — 8, corresponding to G given 
by HI L)| l and the illustration in Fig. [DJ The parameters of 
the idealized system are again a = b = 0.2, and c = 7, 
with a coupling strength of g — 0.79. It is remarkable 
that despite the small number of nodes, the actual lo- 
calized eigenvector and eigenvalue agree well with i|13|) 
(A = 5.334 ... and^ = -0.334 . . . ). 
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FIG. 9: Arrangement of the dynamical units in a ring with 
the strength of the connection between nodes 4 and 5 doubled. 
The matrix G corresponding to this network is in Eq. 1101 . 
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FIG. 10: x$ -14 asa function of time for N = 8 Rossler oscil- 
lators in a ring with the strength of the connection between 
nodes 4 and 5 doubled. The coupling strength is g = 0.79. A 
desynchronization burst starts approximately at t m 9000. 



In Fig. ^]]we show x$ — X4 as a function of time. As 
in the short wavelength case, the burst is not very sharp 
due to the small magnitude of the transversal Lyapunov 
exponents. Nevertheless, it can be seen that the differ- 
ence X5 — X4 increases approximately at t — 9000 and 
returns to a relatively small value after reaching values 
considerably above the average. 

In Fig. Illb we plot the difference between the x co- 
ordinate of node j and its mean over all nodes, Xj — x, 
where x = J2j= 

function of the oscillator in- 
dex j, for t = 8750 (open triangles), t — 9000 (open 



8 



circles), and t = 9250 (open squares). In Fig. II lb we 
show the localized eigenvector of the Laplacian G found 
numerically. As discussed before, the desynchronization 



Xj ~ X 
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FIG. 11: a) Xj — He for t = 8750 (open triangles), t = 9000 
(open circles) , and t — 9250 (open squares) , for the configura- 
tion in Fig. [§] The burst develops with the spatial pattern of 
the localized eigenvector in Fig. lllb . b) Localized eigenvector 
of matrix G in Eq. H1O0 . 



burst affects mainly nodes 4 and 5 (those which share the 
strengthened connection) and the ones adjacent to them. 
Nodes 1,2,7 and 8, however, maintain approximate syn- 
chronization during the burst. 

In Fig. ^] we show the mode weights corresponding to 
the x coordinate as a function of time. The top curve cor- 
responds to [774]! (for the localized mode), and the curves 
close to the horizontal axis to [r]k] x , k ^ 4, for the other 
modes. (The degeneracy of the eigenvalues is broken by 
the strengthened connection, so we do not combine [rik]x 
and [r)N-k]x as before.) Confirming the qualitative sim- 
ilarity between the eigenvector and the spatial pattern 
of the desynchronization burst observed in Fig. ^2 the 
weight corresponding to the localized eigenvector is seen 
to be dominant during the period of time in which the 
burst occurs. 




8500 9500 1 0500 

t 

FIG. 12: [fjk]x as a function of time for k — 4 (top curve) 
corresponding to the localized mode, and for k 7^ 4 (bottom 
curves, close to zero), corresponding to other modes. In the 
burst, the localized mode is excited first and only after some 
time are the other modes also somewhat excited. The local- 
ized mode is dominant during the burst. 



IV. EFFECTS OF THE MISMATCH SPATIAL 
PATTERNS 

In this section we will discuss the effects that the mis- 
match spatial patterns have on the development of the 
desynchronization bursts. For these purposes, it will be 
convenient to rewrite Eq. in the form 

N 

Xi = FiX,) - nY^C.jr. Xy: + Qi(Xi), (14) 
i=i 

where F(Xi) = F(X i: ~p) with ~p = ^2j=iMij and 
Qi(Xi) — F(Xi,/j,i) — F(Xi). The term Qi represents 
the effect of the mismatch and is assumed to be small. 
As before, we linearize around the synchronous state to 
get 

N 

h = DF{ S )e, - 9 J2GijDH(s)e j + Q t {s), (15) 
3=1 

where we have discarded terms of order Qe. With the 
previous notation and Q — [Qi, Q2, ■ ■ ■ Qn], we obtain 
after the diagonalization 

f) k = (DF(s) - g\ k DH{sj) m + (QL) k , (16) 

where (QL) k is the fc'th column of the dx N matrix QL. 
In the ring with equal coupling along each link, the diago- 
nalization procedure corresponds to a Fourier transform. 
In this case, we see that the mismatch affects the different 
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modes according to the weigth, (QL)k, of this particu- 
lar mode in its Fourier expansion. In other cases, for 
example in the localized eigenvector, the strength of the 
mismatch affecting the localized mode is proportional to 
the weigth of the localized eigenvector in the eigenvector 
decomposition of the mismatch. We will now discuss two 
applications of these results. 



A. Amplification of mismatch patterns when 
modes with the same eigenvalue burst 

We have shown that the modes of the mismatch force 
the corresponding modes of the deviations from the syn- 
chronous state. When bubbling induced bursting is ex- 
pected, the size of the mismatch determines the average 
time between bursts [2lJ. Thus, the size of the mismatch 
component in mode k determines the average interburst 
time when that mode is in the bubbling regime. 

When the spectrum of the matrix G is degenerate, the 
spatial modes of the mismatch play an extra role. All 
the modes sharing the same eigenvalue A have the same 
stability properties, and thus, when the corresponding 
value gX is in the bubbling zone, all eigenvectors with 
this eigenvalue are equally likely to appear. The only 
difference between these modes is the strength with which 
they are forced, which is determined by the mismatch 
component in that mode as shown in Eq. (|16|l (or, if 
noise is present, by the noise component in that mode). 

An example of this situation is the ring with connec- 
tions of equal strength in the long wavelength bursting 
scenario. Since the ring is invariant with respect to rota- 
tions, the phase of the long wavelength oscillations can 
not be determined only from the network and dynam- 
ics part of the problem. The two modes with the longest 
wavelength (corresponding to sinusoidal and cosinusoidal 
oscillations) have the same eigenvalue. It is the mismatch 
that in this case determines the phase of the long wave- 
length burst. 

We will show how one can determine the phase of the 
long wavelength desynchronization burst in the case of 
coupled Rossler systems in a ring with equal coupling 
along each link. For this system, the mismatch vector 
Qj(Xj) is given by 




Here 5a = [60,1,62, ■ ■ ■ ,$n] an d similarly for 5b, 5c, and 
y, z are the trajectories around which the linearization 
was made. 

We consider the case in which mismatch in one pa- 
rameter is dominant, for example a. The mismatch in 
the parameters b and c will be assumed negligible com- 
pared with that in a, so that 5b, 5c <C 5a. In this case, 
only the second component of (|18|l is of relevance. Thus 
modes rji and ?7jv-i are excited with a strength propor- 
tional, respectively, to ^(Sa) and !FN-i{5a); sec <|16[1 . 
The magnitude of r/k will be proportional to Tk{5a), and 
thus the excitation of the long wavelength mode (which 
is the only one for which perturbations grow) is propor- 
tional to 



.Myosin ( ^ 



oc sin 



:F\ _i (<VM cos ( '^j- 



2irj 



N 



(19) 



(20) 



(17) 



where tan</> = Tn -\{5a) / J-\(5a). 

We now show results of numerical simulations illustrat- 
ing the above. The parameters N and g will be as in the 
long wavelength example in the previous section. We use 
the same random set of perturbations used in that exam- 
ple. As described above, we obtained the phase <f> of the 
long wavelength component of the vector 5ai. In Fig. 1131 
we plot yj — Uj-i for different times during a burst (filled 
symbols). In the same Figure, we plot a scaled version 

of sin (2^ + 0)— sin ^ 27r( -j' 2 ~ 1 - ) + (j^j (open circles). The 

phase of the desynchronization burst is in agreement with 
that of the long wavelength component of the mismatch. 

When the mismatch affects predominantly one param- 
eter as in this case, the phase of the bursts can be pre- 
dicted as described above. When mismatch in different 
parameters is comparable, the phases of the long wave- 
length modes of the different parameter mismatches com- 
pete and the bursts develop with one of these phases or 
with a combination of them. 

It must be emphasized that this analysis is possible 
only when there is a degeneracy of the eigenvalues. For 
example, the location of the localized burst can not be 
determined in this way, as it is fixed in the position of the 
strengthened link. In this case, the mismatch component 
in the localized mode would only affect the average time 
between bursts. 



a und similarly for 



bj and 5c j. 



We 



where is the normalized 



where 5aj — aj — 
define T k {u) = Y,J=i 

j'th component of the k eigenvector described at the be- 
ginning of Section III. With this convention, the term 
(QL)k in equation i|ltj|) is given by 





(QL) k = ( yFk{5a) 

T k {5b) - zT k (5c) 



(18) 



B. Artificial supression of unstable modes using 
knowledge of the mismatch 

We will now discuss another consequence of Eq. Ultifl . 
We imagine a situation where we are given a number of 
nearly identical oscillators that we are to connect in a 
network which we desire to be in synchronism as much 
as possible. Furthermore, we imagine that, through mea- 
surements made individually on each oscillator, we are 
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y>- yj-i 




FIG. 13: ?/j — Dj-i for different times during a burst 
(filled symbols), and a scaled version of sm(^j-+(f\ — 

sin(^=^- + A with <t> as given in the text (open circles). 
The phase of the burst spatial pattern coincides with the 
phase of the long wavelength component of the mismatch. 



plot X\ — X2 as a function of time for configuration ai arge 
(top curve) and for configuration a sma ii (bottom curve). 
The difference x\ — x^ is much smaller in the former case 

X ] ~ Xq 



'fi 




a 

LA 


arge 


Aw 






_1 1 








I 


"mall 





5000 10000 

t 



FIG. 14: si -i2 as a function of time for a configuration of 
oscillators with a large (top curve) and with a small (curve 
closer to zero) short wavelength component of the mismatch. 
The quality of the synchronization is much better in the sec- 
ond case. 



aware of the amount of mismatch in each oscillator. The 
question we address is this: Using our knowledge of the 
individual mismatches, how should we arrange the os- 
cillators in the network so as to best supress the fre- 
quency of desynchronism bursts? To answer this ques- 
tion, we note that, according to the previous discussion, 
we should reduce the mismatch component in the mode 
which is in the bubbling region. Since the size of the mis- 
match affects the average interburst time [2l]|, reducing 
this component is desirable if one wants to improve the 
quality of the synchronization. This can be accomplished 
by judiciously arranging the dynamical units so that the 
fc'th mode of the mismatch is minimized when the corre- 
sponding value g\k is in the bubbling region. For exam- 
ple, to supress long wavelength bursts, one may arrange 
the units so that the parameter errors alternate above 
and below the mean. To supress the localized bursting 
described in the previous section, one could arrange the 
units so that those with the more similar parameters are 
the ones in the region of the strengthened connection. 

As a concrete example, we test this idea using simula- 
tions for the case of short wavelength bursting presented 
in the previous Section. We again assume for simplicity 
that mismatch in the parameter a is dominant. We gen- 
erate random perturbations in the parameter a within 
a ±6% range of the value a — 0.2, as explained in the 
previous section. With this set of parameters given, we 
set up the dynamical units in the ring using two different 
permutations of their positions. One of them (a s ) has a 
smaller and the other (a; ) a larger short wavelength com- 
ponent Tiia) than the original random sequence. The 
ratio T^ai) / Ti(a s ) is approximately 15. In Fig 1141 we 



than in the latter, roughly by a factor of 15, as can be ex- 
pected from the ratio T^{ai arge ) / '^(a small) ■ This quali- 
tative example illustrates how one can use knowledge of 
the mismatch to supress undcsircd instabilities. 

V. SPATIAL PATTERNS OF DEVIATIONS 
FROM THE STABLE SYNCHRONOUS STATE 

So far, we have concentrated in the case in which the 
value of gXk is in the bubbling regime for one mode k and 
in the stable regime for the other modes, so that desyn- 
chronization bursts occur sporadically. As we have seen, 
these bursts present spatial patterns on the network. 

If synchronization is desired, one would might try to 
avoid the bubbling regime by designing the network and 
adjusting the coupling strength so that all the modes lie 
in the stable zone. One would also strive to reduce the 
mismatch, but as mentioned before, there are practical 
limitations on how much one can make the oscillators 
exactly the same. 

If ^(gXk) is negative for all modes (indicating transver- 
sal stability of the synchronous state) one can have, de- 
pending on the degree of transversal stability, fair syn- 
chronization even with relatively large amounts of mis- 
match. If one is to operate under such conditions, it is 
important to know the characteristics of the deviations 
from the synchronous state. 

Thus we ask in this scenario: How large are the spatial 
patterns of the deviations from the synchronous state, 
and how does this depend on the mismatch and on the 
degree of transversal stability? 



11 



The spatial modes of these deviations obey Eq. (|16fl . 
In the absence of the term (QL)k, the zero solution is 
stable, and typical perturbations from it decay, having 
a negative Lyapunov exponent given by hk = ^(gXk)- 
The first term in the right hand side of Eq. (|16fl can 
be thougt of as a damping term with a damping rate 
given by hk, and the second term, {QL)k, as a forcing 
term. Since we are considering the stable case, these 
two factors, on average, cancel each other. By definition, 
the Lyapunov exponent for the system without mismatch 

is given by hk = ( Vk ^ DF ^ x ^ DH ^ !rik ^ where the angle 
brackets indicate time average. Assuming a solution rjk 
of the system with mismatch to yield the same value of 
this time average, we left multiply Eq. (|16[) by rj^ |%| -2 
and average to obtain 

M « <^^> ~ <^j4 (21) 
1*7*1 W 

where the angle brackets indicate time average. This 
leads to the following rough estimate, 

,v (\m) k \) r99 , 

\h,k\ 

As an example we consider Rossler units in a ring with 
all connections of equal strength. We choose N = 8, 
g = 0.6 (^(g\k) < for all values of k). Furthermore, 
we add a random perturbation to the parameter a of each 
oscillator chosen uniformly from within a ±0.1% range of 
a = 0.2. 

In Fig. E| we show, for k = 1,...,7, the quanti- 
ties (|i7fc|) (squares), (\(QL) k \) (triangles), and 
(stars) . The magnitudes of the forcing term for the differ- 

□ :<\>, A\<(QL) k >, 0:<(QL) fc > 




1 2 3 4 5 6 7 

k 

FIG. 15: {\gu\) (open squares), (\(QL)k\) (open triangles), 
and <l(Q . £)fcl> (open circles) for N = 8, g = 0.6, ft = 1, . . . , 7. 
The forcing term (open triangles) roughly determines the re- 
sponse (open squares). The corrected forcing term (open cir- 
cles) matches well the response (open squares). 



ent modes ((\(QL)k\)) span roughly two orders of magni- 
tude, and the magnitude of the response {(\g k \)) looks 
roughly proportional to the latter. When the forcing 
term is corrected by dividing it by the magnitude of the 
corresponding Lyapunov vector \hk\, the resulting quan- 
tity ( <l( f^] fcl> ) matches very well the observed response. 



VI. CONCLUSIONS 

We have studied the stability properties of the syn- 
chronized state in a network of coupled chaotic dynami- 
cal units when these have a small heterogeneity. We have 
shown that when the dynamical units that are coupled 
in a network are sligthly different, the synchronized state 
can be interrupted by large infrequent desynchronization 
bursts for some values of the parameters. The range of 
the parameters for which this phenomenon is expected 
can be obtained by performing a master stability func- 
tion analysis of the chaotic attractor and of the periodic 
orbits embedded in it. 

The desynchronization bursts are induced by the bub- 
bling phenomenon, and have spatial patterns on the net- 
work. These spatial patterns can be predicted from the 
eigenvectors of the Laplacian matrix G and the master 
stability functions mentioned above. We showed exam- 
ples illustrating the development of bursts with spatial 
patterns. One of our examples showed that the strength- 
ening of a single connection might destabilize the nodes 
near this connection, while leaving the rest of the network 
approximately synchronized. 

Direct measurement of the parameter mismatch in the 
elements of a network might prove useful. We discussed 
how this knowledge could be used to reduce the frequency 
of bursts and to predict the relative weights of different 
spatial patterns in a burst. We also discussed how one 
could, from knowledge of the mismatch and of the mas- 
ter stability function, describe the spatial patterns and 
magnitude of the deviations from the synchronized state 
when the synchronization of the corresponding identical 
unit system is robust. 

We emphasize that although we did not discuss the 
effects of noise, the phenomenon described in this paper 
also occurs for noisy identical oscillators. Desynchroniza- 
tion bursts with spatial patterns are expected for noisy, 
identical oscillators if one has them for noiseless, non- 
identical oscillators. The difference is that the param- 
eter mismatch is always 'frozen', in the sense that the 
mismatch is always the same for each oscillator, whereas 
for noise this is not the case. 
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